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ABSTRACT 

The magnetic 0-star HD 191612 exhibits strongly variable, cyclic Balmer line 
emission on a 538-day period. We show here that its variable Ha emission can be well 
reproduced by the rotational phase variation of synthetic spectra computed directly 
from full radiation magneto-hydrodynamical simulations of a magnetically confined 
wind. In slow rotators such as HD 191612, wind material on closed magnetic field 
loops falls back to the star, but the transient suspension of material within the loops 
leads to a statistically overdense, low velocity region around the magnetic equator, 
causing the spectral variations. We contrast such "dynamical magnetospheres" (DMs) 
with the more steady-state "centrifugal magnetospheres" of stars with rapid rotation, 
and discuss the prospects of using this DM paradigm to explain periodic line emission 
from also other non-rapidly rotating magnetic massive stars. 
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1 INTRODUCTION 



Shortly after Donati et al. ( 2006 1 detected a strong magnetic 



field i n the Galactic Of?p star HD 191612, [Howarth et al.| 
(|2007[) demonstrated that the variable equivalent widths 



of its optical Balmer and He I lines (e.g., Walborn et al 



2003) can be accurately phased according to a 538-day pe- 



riod, where in particular the outstanding Ha variation shows 
strict periodicity. Since this period is unrelated to the much 
longer orbital period Porb = 1542 d of HD 191612 and its 
binary companion ( Howarth et al.|[2007 l, rotational modu- 
lation of a magnetically confined wind seems the most likely 
origin for the variability, as already suggested by |Donati| 



et al. ( 2006 1. But in contrast to centrifugally supported mag- 



netosphere models, which have been successfully applied to 
Balmer line variability in rapid rotators such as the B star 
a Ori E ( Townsend et al!]|2005 l, it is not clear how a very 
slow rotator such as HD 191612 can sustain a magnetosphere 
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with sufficient accumulation of wind plasma to explain the 
strong and periodic Balmer emission. 

To reproduce the Ha variation of HD 191612, [Howarth] 



et al. (20071 suggested two geometrical toy models. One of 



these was indeed inspired by the plasma distribution qual- 
itatively expected from a magnetically confined wind; it is 
a tilted, limb-darkened, geometrically thin disc, where the 
sum of observer inclination i and obliquity /3 (the angle be- 
tween the rotation and magnetic axes) must be i -|- /3 ~ 100° 
for the Ha modulation to be fit. 



Wade et al. (20111 recently analysed Stokes V spectra 



of HD 191612. Assuming a dipole oblique rotator, these au- 
thors derived i -f- /? = 95 ± 10°, and by matching electron 
scattering modelling to the observed photometric variabil- 
ity further obtained i ^ 30° . A tentative reference geometry 
i = 30° and /3 = 67±5° was then suggested from speculating 
that the orbital and spin angular momenta of HD 191612 be 
aligned, and a surface dipole (polar) field B^ = 2450 ±400 G 
derived. 

This Letter examines to what extent full radiation 
magneto-hydrodynamical (MHD) simulations of a magnet- 
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Observations 

Spherically symmetric models 1 




Velocity [km/s] 



Figure 1 . Observed jHowarth et al.|2007| and synthetic Ha spec- 
tra during piiases close to minimum and maximum. Synthetic 
FASTWIND spectra are computed for two different mass-loss rates 
under the assumption of spherical symmetry (see Sect.[2]l. 



Table 1. Summary of stellar, wind, and magnetospheric param- 
eters of HD 191612 ( [Howarth et al.|2007|[Wade et al.|20lT| |. 



Name 



Parameter Value 



Effective temperature 


TLff 


35 000 ±1 000 K 


Surface gravity 


logg 


3.5 ±0.1 


Stellar radius 


R. 


14.5 i?0 


Helium abundance 


n-no/n-a 


0.1 


Terminal speed 


Vao 


2700 km/s 


Mass-loss rate and 


M^i 


1.6 X lO-'^Me/yr 


clumping factor 






Surface polar mag- 


Bi 


2450 ± 400 G 


netic field 






Obliquity and ob- 


H + i 


95 ± 10° 


server inclination 







ically confined wind, along with detailed radiative transfer 
calculations, can actually reproduce HD 191612's observed 
Hq variability, under the wind, magnetic, and geometric 



constraints derived by Howarth et al. ( 2007 1 and Wade et al, 

mn\. 



2 Hq IN A SPHERICALLY SYMMETRIC 
WIND MODEL 

To set the stage, we first compute synthetic Ha profiles for 
two different mass-loss rates using the spherically symmet- 
ric, unified (photosphere+wind) NLTE (=Non Local Ther- 
modynamic Equilibrium) model atmosphere code fastwind 
( Puis et al.|2005 l, taking stellar and wind parameters from 
[Howarth etal.| ( |2007| ) (Table [T|. Fig. [l] confronts such mod- 
els with the observed Ha line profiles during minimum and 
maximum phases. 

The model-fit during minimum phase is acceptable. In- 
deed, it is equivalent to that of |Howarth et al 
from which M 



/hi 



(20071, 



1.6 X 10 M0/yr was derived for 
HD 191612. Here /d = j^ ^ 1 is the wind clumping fac- 
tor, with average mass density (p). /d enters the analysis 
because Hq is a recombination based line in O stars, which 
means that its line strength is greater in a small-scale struc- 
tured ('clumped') wind than in a smooth wind with the same 
mass- loss rate (e.g., Sundqvist et al.|2011|. 



In an attempt to also match the maximum phase, we 
next calculated a test-model with 5 times higher mass-loss 
rate. However, Fig. [T] clearly shows that this synthetic Ha 
profile is much broader than the observed one. This indicates 
that the variable Hq emission does not stem from a variable 
global mass-loss rate accompanied by a spherically symmet- 
ric velocity field such as the v — iioo(l — ^/r) field assumed 
in FASTWIND, here with /3 = 1. Rather the variability might 
be caused by a confined region of wind material with high 
density and low velocity; such confinement may indeed stem 
from the strong magnetic field in HD 191612 channeling its 
radiatively driven wind outfiow to form a stellar magneto- 
sphere. We now describe our efforts to model this hypothe- 
sized structure. 



3 SIMULATIONS 



3.1 Modelling a dynamical magnetosphere 



Following the general procedure outlined by |ud-Doula fc| 
|Owocki| ( |2002[ ), we compute a 2-D radiation MHD wind 
simulation of HD 191612, assuming a dipole magnetic field. 
Hydrodynamical variables are specified on a standard, right- 
handed spherical grid (r, 9,(j)), defined relative to a Cartesian 
set {x,y,z), where we assume symmetry in <f). The energy 



equation is treated as by Gagne et al. ( 2005 I and the ra- 



diation line force is calculated within the Sobolev approx- 



imation using standard CAK ( Castor et al. 1975 1 theory. 



Since the rotation of HD 191612 is extremely slow, the in- 
ferred period of 537.6 days ( [Howarth et al.|2007[ ) implies an 
equatorial rotation speed «rot ~ 1.4 km/s, we may neglect 
rotational effects on the dynamics (and thus use the same 
simulation for any choice of obliquity /3). 

The effectiveness of the magnetic field in channeling the 
stellar wind outflow may be characterized by the ratio of 
magnetic to wind kinetic energy density. 



B^ 



pv^/2 



TT (r/R*)~ 



(1) 



where the second equality deflnes the so -called 'wind confine- 
ment parameter' r;* = BlRl/{Mvrx) (ud-Doula & Owocki 



20021, with _B* the dipole equatorial surface field strength. 
> 1, the dipole Alfven radius Ra « ?)* i?, , at which 



If??* 

the magnetic and wind energy densities are equal, is located 
away from the stellar surface, allowing then for some wind 
material to be channeled along closed loops towards the 
magnetic equator (see Fig. pi. But the much steeper ra- 
dial decline of the dipole magnetic energy density (~ 1/?"^) 
than the wind kinetic energy density (~ l/r^), means that 
at large enough radii the wind will always force the field 
lines to open up and essentially follow the radial wind flow. 
The simulation here assumes strong conflnement, 77* — 
50, in accordance with the magnetic field strength recently 
derived by [Wade et al.[(2011[) and the wind parameters de- 
rived by Howarth et al. (20071, adopting /d = 1. It is well 
established that the winds of hot, massive O stars are indeed 
clumped (see [Sundqvist et al.|201l"| for a recent review). But 
a theoretical development of such stochastic, small-scale in- 
homogeneities, as caused by the strong instability inherent 
to line-driven winds (e.g., [Owocki et al.|[T988l ), requires a 
non-Sobolev treatment of the radiation line force, and has 
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yet to be implemented within any MHD simulation. How- 
ever, in terms of the Ha modelling that is the focus of this 
paper, we are still effectively modelling M^fjik (see Sect.|2|, 
but simply neglecting any dynamical effects such stochastic, 
small-scale structures might have upon the large-scale wind 
structure imposed by the magnetic field. 

The upper panels of Fig. [2] plot the density squared 
of two simulation snapshots. They illustrate how below 
r ~ Ra ~ 2.7-R*, the wind does indeed become trapped 
by the closed field-line loops, whereby the material is pulled 
back by gravity onto the star over a dynamical time-scale. 
But a key point here is that, despite the very dynamical 
behaviour, the transient suspension of material within such 
closed loops still results in a wind region, in the vicinity of 
the magnetic equator, that statistically is overdense (Fig. [2| 
lower-left panel). Further, as a result of the colliding wind 
material at individual loop tops, this overdense region is 
also characterized by very low velocities (Fig. [2] lower-right 
panel) , in qualitative agreement with the narrowness of the 
observed Ha emission discussed in Sect. 121 

The structures predicted by these simulations are phys- 
ically distinct from those predicted for rapidly rotating mag- 
netic stars with R\ > J?k ([Townsend fc Owocki| |2005{ 



ud-Doula et al. 2008), where _Rk 



Townsend et al.||2007 

iJhotpUcrit)^n* is the Kepler co-rotation radius for critical 
rotation speed Wcrit- For such stars, the centrifugal forces can 
support any trapped material above Rk, allowing then the 
magnetically confined wind to accumulate material and form 
a centrifugal magnetosphere (CM). In contrast, the char- 
acteristic structure described above, appropriate for slowly 
rotating massive stars with Rk > -Ra > R*, instead estab- 
lishes the concept of a dynamical magnetosphere (DM) (see 
also [Petit et al.|2011 1 . 



In hot coronae from the sun and some magnetically ac- 
tive cool stars, there are analogous examples of regions of 
dynamical infall ("coronal rain"; e.g. Elbe et al. 19991 or 



centrifugally supported prominences ([Collier Cameron et~al 



2003 Jardine & van Ballegooijen 20051, fed largely by the 



transient eruptive propulsion of stellar flares. By contrast, 
hot-star magnetospheres are fed by the quasi-steady wind 
upfiow driven by the star's radiation, allowing for persistent 
Balmer emission that has been monitored over multi-year 
time-scales spanning many rotation periods. 

Note that even rapid rotators will have a DM compo- 
nent at r < Rk- But the Ha emission contribution from this 
part will be insignificant because of the much higher densi- 
ties at Rk < r < Ra (see, e.g.. Fig. 7 in Townsend et al. 
[2007[ ). These higher densities stem from the much longer 
accumulation time-scale associated with a CM (typically 
months/years, see Appendix A in Townsend & Owocki 2005 1 
than with a DM (typically hours, the dynamical time-scale). 
But a star such as HD 191612, with Ra « 2.7i?* << -Rk ~ 
557?*, only has a DM contributing to the Ha emission. So 
whereas rapidly rotating magnetic B-stars can indeed show 
substantial Balmer line emission, as observed in e.g. a Ori E, 
the short accumulation time-scale of a DM requires the rela- 
tively high mass-loss rate of an O-star to produce observable 
Ha emission in these slowly rotating stars. 




Figure 2. Contours of the density-squared for two different snap- 
shots of the MHD wind simulation (upper panels), and of time- 
averaged density-squared (lower left) and radial velocity (lower 
right). Time averages are calculated from > 100 snapshots taken 
well after the simulation's initial state. The Alfven radius is here 
located at r Ri 2.7 R^, whereas the Kepler co-rotation radius is at 
r Ri 55-R*, i.e. outside the range of the plots. 



3.2 Radiative transfer 

To model the observed variation of Balmer emission, we 
compute synthetic Ha flux profiles directly from the MHD 
simulations by solving the formal integral of radiative trans- 
fer in a 3-D cylindrical coordinate system {p,^,z'). This 
system is aligned toward the observer by rotating the stel- 
lar system (r, 0, (j>) by an angle a about its y-axis, so that 
cos a = z-z . The angle a thus defines the observer's viewing 
angle with respect to the magnetic pole. For given /3 and i, 
we have 



cos a — sin /3 cos $ sin i + cos /3 cos i, 



(2) 



which then readily gives the observer's viewing angle as func- 
tion of rotation phase <I>. We note that even though the MHD 
models are 2-D, the radiative transfer must be performed in 
3-D, as the axial symmetry is broken for any observer with 
I cosa| 7^ 1. 

The transfer equation is solved only in the wind, with 
a pre-specified photospheric profile P^ as a lower boundary 
condition, taken from NLTE model atmosphere calculations 
assuming negligible wind contamination. While not truly 
self-consistent, this procedure has been shown to be very ac- 



curate for Ha line profile calculations in 1-D smooth (Puis 



et al. 11996 1 as well as multi-dimensional clumped ( Sundqvist 



et al. 12011 1 O-star wind models without magnetic fields. 



The monochromatic optical depth along a ray is 



(3) 



where Xi^ is the frequency dependent opacity per unit length. 
The opacities are calculated assuming an optically thin con- 
tinuum and occupation numbers for the Ha atomic lev- 
els i given in terms of the NLTE departure coefficients 
fei ~ Ui/rit. Here n* is the occupation number of level i in 
LTE with respect to the ground state of the next ionization 
state (e.g., [M"ihalas[[T978[ ). The line profile is a Gaussian of 
Doppler width A/^d, set by the local wind electron temper- 
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ature Tc, and centered at zero co-moving frame frequency 

3::cmf = a^obs — z' ■ u/uoo, where x = {v/vq — l)c/voc- 
The emergent intensity for a given ray then is 



/,. = Puloe 



+ 



Sv{Ti,)e ^"dr^, 



(4) 



where /q is the stellar photospheric continuum intensity, Sv 
the NLTE line source function, and r^ the optical depth 
integrated over the complete ray. /o is taken from fastwind 
model atmospheres for rays that intersect the stellar core 
and set to zero otherwise. Su is fixed by the Ha departure 
coefficients and the local wind electron temperature. The 
emergent flux is then, finally, obtained by integrating the 
emergent intensity over the projected stellar disc. 
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Figure 3. Synthetic Ha profiles for 8 randomly selected snap- 
shots during our simulation run (dashed lines), and a mean profile 
calculated from them (solid line). 



3.2.1 Electron temperatures and Ha occupation numbers 

As described, the Hq synthesis problem requires estimates 
of Te and the hydrogen departure coefficients. But the en- 
ergy equation as treated in the MHD simulations described 
in Sect. |3.1| yields only a rough approximation of the wind 
temperature balance, with the local temperature artificially 
never allowed to drop below a certain floor value (on the 
order of the stellar effective temperature). In our Ha cal- 
culations, we therefore estimate the wind temperature bal- 
ance using the results of a spherically symmetric fastwind 
model, except in regions shock-heated to T^ > lO^K, where 
we set the Ha opacity and source function to zero. 

To consistently calculate NLTE departure coefficients 
for a full multi-D MHD wind simulation is a daunting task, 
well beyond the scope of the present paper. However, for 
now we take advantage of the fact that hydrogen is almost 
fully ionized in typical O star winds. The Ha line forma- 
tion is then controlled by recombination, a thermal process, 
and the participating atomic levels are therefore very close 
to LTE with respect to ionized hydrogen. For HD 191612, 
FASTWIND calculations show deviations smaller than a fac- 
tor of two, a typical number for most 'normal' O-star winds 
without strong magnetic fields (Puis et al.||1996 Sundqvist 



et al. 2011 1. As a first approximation then, we here take 



the simplest approach possible and assume LTE conditions, 
whereby b^ — 1. 



4 Ha VARIABILITY IN A DYNAMICAL 
MAGNETOSPHERE 

We calculate an Ha line profile for more than 100 snapshots 
of the 2-D MHD simulation of HD 191612. Fig. |3]shows that 
such profiles are highly variable, mimicking the wind's dy- 
namical nature. To obtain a mean profile at each phase, 
we average over such profiles for ~100 randomly selected 
snapshots; this is intended to be a simple proxy for the 
real 3-D nature of the wind dynamics, effectively using the 
complex and non-linear variations in time in our 2-D sim- 
ulation to mimic the expected variations in azimuth in a 
full 3-D modelr] While approximate, this seems a reason- 
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Figure 4. Observed (ephemeris from [Howarth et al. | |2007[ l 
and synthetic Ha dynamic spectra, as functions of rotation 
phase and shown over two cycles. Synthetic spectra assume 
P = i = 50° and are calculated as described in text. The 
model spectra in the lower panel have further been convolved 
with an isotropic 'macro-turbulence' of 100 km/s. 



For this type of line transfer, which often is optically thick, 
such post-averaging of the line profiles computed for many time 
snapshots is more realistic than the simple pre-processed, time- 
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Figure 5. Upper panel: Observed (filled circles) and simulated 
Hcf equivalent widths as functions of rotation phase, for the three 
geometries j = ^ = 50° (red, solid), i = 30°, /3 = 70° (dotted), 
and i = 10°, /3 = 90° (dashed). The blue filled circles represent 
the higher quality data points used in Fig. |4] Lower panels: Pro- 
jected surface area contours of the normalized emergent intensity 
at line center, for observers located along the magnetic equator 
(left) and pole (right). Cartesian coordinate x' on the abscissa 
and y' on the ordinate. 



able first approach to account for the limited lateral coher- 
ence or synchronization of an actual 3-D magnetized wind. 
Also, each phase is constructed by further averaging the two 
phases having equal <l>'s when reversing the magnetic poles, 
to ensure the expected long-term north/south symmetry of 
our simulation. The procedures described above effectively 
smooth out most of the short-time variability of our simu- 



lation, in agreement with the observations (Howarth et al 



2007). 



Fig. |4] compares observed and synthetic time-averaged 
dynamic Ha spectra, plotted as functions of rotational phase 
assuming P = i — 50° . This is consistent with the P + i — 



95 ± 10° derived by |Wade et al.| ( |2011| ), but differs slightly 
from the i = 30° adopted there (see further below). The 
observed general trends, with peak flux at phase and an 
extended minimum around phase 0.5, are both well repro- 
duced. The flux variations are caused by differences in the 
projected surface area of overdense Ha emitting material as 
the observer changes viewing angle when the star rotates. 
Fig. [5] demonstrates this by plotting the Ha emergent inten- 
sity (surface brightness) at line center for observers located 
along the axes of magnetic pole and equator. The flgure 
clearly shows how the flux, which is just the integral of the 
intensity over this projected area, is much higher for the 
observer along the polar axis. 

The large observed Ha variability puts rather tight con- 
straints on the system's geometry. The equivalent width 
curves in Fig. |5] directly refute very low values of i, but 
also show that the i = 30°, /3 = 70° assigned as a tentative 



reference geometry by Wade et al. (20111 results in some- 
what weaker variation than the fS — i = 50° adopted here. 



This is simply because, for a given sum i -\- P = 100°, an 
observer a,t i — 30° never looks closer to the magnetic pole 
than a ~ 40°, whereas for an observer sd, i — 50°, a spans 
the entire range from pole to equator, and back again, in 
one rotation period, thus resulting in larger flux variations. 
There are discrepancies, of course. Whilst the observed 
and simulated equivalent width curves qualitatively agree 
well, the simulated variation is quantitatively somewhat too 
low. These deviations could however be remedied if the DM 
were more concentrated, which would result in a larger sur- 
face brightness difference between pole and equator (Fig.[5|. 
Such stronger wind confinement could occur from either a 
lower mass-loss rate or a stronger magnetic field, where the 



former choice seems more likely (due to clumping. Sect. 3.1 1, 
since the i — 50° adopted here actually would result in 
a slightly reduced magnetic field strength as compared to 



that derived by Wade et al. (20111, who adopted i — 30 



Another possibility is of course that these slight discrepan- 
cies simply are related to insufficient assumptions for the 
wind electron temperature structure and/or the hydrogen 



occupation numbers (Sect. 3.2.11 



In addition, the velocity dispersion in the models is too 
low, predicting narrower and sharper peaked profiles than 
observed (Fig. 14] upper panel). To illustrate this further, 
the lower panel of Fig. H] displays the same line profiles as 
the upper panel, but now with the model profiles convolved 
with an isotropic Gaussian 'macro-turbulence' of lOOkm/s. 
Inspection shows that the new profile shapes are significantly 
improved. We suspect that this required macro-turbulence 
is an artefact of missing wind dynamics, since, generally, the 
added degree of freedom in a 3-D simulation should result in 
larger velocity dispersion than in the 2-D models employed 
here. Though very challenging, 3-D MHD models are cur- 
rently being developed by one of us (A. ud-Doula), and will 
be reported in a future paper. Such 3-D simulations will then 
also quantitatively test our simple approach here of averag- 
ing radiative transfer results for 2-D simulation snapshots 
to approximate the real 3-D dynamical wind structure. 

Finally, notwithstanding the foregoing comments, there 
is surprisingly good overall agreement between the mod- 
els and observations, which strongly supports that the DM 
model captures the key physics responsible for the Ha vari- 
ability. 



5 DISCUSSION AND CONCLUSIONS 

We have demonstrated that radiation magneto- 
hydrodynamical simulations of a confined wind, together 
with detailed radiative transfer modelling, reproduce well 
the distinct periodic Ha emission observed in the magnetic 
O-star HD 191612. We interpret this within the context of 
a dynamical magneto sphere (DM), wherein the rotationally 
modulated spectral variations are results of a statistically 
overdense, low velocity wind region around the magnetic 
equator. 

While applied here only to HD 191612, the DM model 
may also well describe optical Balmer line variability in other 
magnetic O stars with 7?k > .Ra, such as 9^ Ori C and 
HD 148937 (IWade et al.||2006l 120111). Indeed, the narrow 



averaged wind density used by [Wade et al.| | |2011[ l to model po- 
larized electron scattering, which is more nearly optically thin. 



Ha emission observed in HD 148937 suggests a line forma- 
tion scenario in a DM, with the only significant difference 
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to HD 191612 then being the stellar and/or magnetic geom- 



etry ( Wade et al.|2011 1, resulting in much smaller spectral 
variations for the former star. In future work, we intend to 
develop further this DM model and apply it to a broader 
sample of magnetic massive stars that show variable Balmer 
emission and are characterized by 7?k > Ra- 
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